import datetime
import os
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
import plotly.graph_objects as go
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
from sklearn.metrics import mean_absolute_error
from prophet import Prophet
import statsmodels.api as sm
from statsmodels.stats.diagnostic import acorr_ljungbox

import warnings
warnings.filterwarnings("ignore")

# Format the date as a string
formatted_date = datetime.date.today().strftime("%B %d, %Y")

text_block = f"""
This notebook was executed on: {formatted_date}.
Any results or analysis reflect the data and conditions as of this date.
"""

print(text_block)

This notebook was executed on: May 07, 2025.
Any results or analysis reflect the data and conditions as of this date.

hourly_df = pd.read_csv("pharma-sales-data/saleshourly.csv")
hourly_df.head()

hourly_df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 50532 entries, 0 to 50531
Data columns (total 13 columns):
 #   Column        Non-Null Count  Dtype  
---  ------        --------------  -----  
 0   datum         50532 non-null  object 
 1   M01AB         50532 non-null  float64
 2   M01AE         50532 non-null  float64
 3   N02BA         50532 non-null  float64
 4   N02BE         50532 non-null  float64
 5   N05B          50532 non-null  float64
 6   N05C          50532 non-null  float64
 7   R03           50532 non-null  float64
 8   R06           50532 non-null  float64
 9   Year          50532 non-null  int64  
 10  Month         50532 non-null  int64  
 11  Hour          50532 non-null  int64  
 12  Weekday Name  50532 non-null  object 
dtypes: float64(8), int64(3), object(2)
memory usage: 5.0+ MB

daily_df = pd.read_csv("pharma-sales-data/salesdaily.csv")
daily_df.head()

weekly_df = pd.read_csv("pharma-sales-data/salesweekly.csv")
weekly_df.head()

# Convert datatype of datum
hourly_df['datum'] = pd.to_datetime(hourly_df['datum'])
daily_df['datum'] = pd.to_datetime(daily_df['datum'])
weekly_df['datum'] = pd.to_datetime(weekly_df['datum'])

hourly_df.set_index('datum', inplace=True)

# Create figure
fig_hourly = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_hourly.add_trace(go.Scatter(
        x=hourly_df.index,
        y=hourly_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_hourly.update_layout(
    title='Interactive Hourly Drug Sales by Category',
    xaxis_title='Datetime (Hour)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_hourly.show()

daily_df.set_index('datum', inplace=True)

# Create figure
fig_daily = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_daily.add_trace(go.Scatter(
        x=daily_df.index,
        y=daily_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_daily.update_layout(
    title='Interactive Daily Drug Sales by Category',
    xaxis_title='Datetime (Day)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_daily.show()

weekly_df.set_index('datum', inplace=True)

# Create figure
fig_weekly = go.Figure()

# Add a line trace for each drug category
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']
for col in categories:
    fig_weekly.add_trace(go.Scatter(
        x=weekly_df.index,
        y=weekly_df[col],
        mode='lines',
        name=col
    ))

# Update layout
fig_weekly.update_layout(
    title='Interactive Weekly Drug Sales by Category',
    xaxis_title='Datetime (Day)',
    yaxis_title='Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1200,
    height=600
)

# Show plot
fig_weekly.show()

# Resample hourly data into monthly total sales per drug
monthly_agg = hourly_df[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']].resample('M').sum()

# Create figure
fig_monthly = go.Figure()

# Add a line for each drug category
for col in monthly_agg.columns:
    fig_monthly.add_trace(go.Scatter(
        x=monthly_agg.index,
        y=monthly_agg[col],
        mode='lines+markers',
        name=col
    ))

# Customize layout
fig_monthly.update_layout(
    title='Interactive Monthly Drug Sales by Category',
    xaxis_title='Month',
    yaxis_title='Total Sales',
    hovermode='x unified',
    template='plotly_white',
    width=1000,
    height=500
)

fig_monthly.show()

# List of drug categories
categories = ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']

# Set weekday order
weekday_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']

# Create subplots for each category
plt.figure(figsize=(16, 16))
for i, cat in enumerate(categories, 1):
    weekday_avg = daily_df.groupby('Weekday Name')[cat].mean().reindex(weekday_order)
    
    plt.subplot(3, 3, i)
    sns.barplot(x=weekday_avg.index, y=weekday_avg.values, palette='viridis')
    plt.title(f'Average {cat} Sales by Day')
    plt.ylabel('Avg Sales')
    plt.xticks(rotation=45)

plt.tight_layout()
plt.show()

# Calculate total average sales for each weekday
daily_df['Total Sales'] = daily_df[categories].sum(axis=1)  # Sum across categories
weekday_avg = daily_df.groupby('Weekday Name')['Total Sales'].mean().reindex(weekday_order)

# Create a single plot for total average sales
plt.figure(figsize=(12, 8))

# Create the bar plot
ax = sns.barplot(x=weekday_avg.index, y=weekday_avg.values, palette='viridis')

# Add labels to each bar
for p in ax.patches:
    height = p.get_height()
    ax.annotate(f'{height:.2f}',
                (p.get_x() + p.get_width() / 2, height),
                ha='center', va='bottom')

plt.title('Average Total Sales by Weekday from All Products')
plt.xlabel('Weekday')
plt.ylabel('Average Total Sales')
plt.xticks(rotation=45)

plt.tight_layout()
plt.show()

# Set plot style
sns.set(style='whitegrid')

# Create subplots
plt.figure(figsize=(16, 12))
for i, cat in enumerate(categories, 1):
    hourly_avg = hourly_df.groupby('Hour')[cat].mean()
    
    plt.subplot(3, 3, i)
    sns.lineplot(x=hourly_avg.index, y=hourly_avg.values, marker='o')
    plt.title(f'Average Hourly Sales for {cat}')
    plt.xlabel('Hour of Day')
    plt.ylabel('Average Sales')
    plt.xticks(range(0, 24, 2))

plt.tight_layout()
plt.show()

# Set plot style
sns.set(style='whitegrid')

# Calculate total average hourly sales
hourly_df['Total Sales'] = hourly_df[categories].sum(axis=1)
hourly_avg = hourly_df.groupby('Hour')['Total Sales'].mean()

# Create a single plot for total average hourly sales
plt.figure(figsize=(10, 6))
sns.lineplot(x=hourly_avg.index, y=hourly_avg.values, marker='o')
plt.title('Average Hourly Sales for All Products')
plt.xlabel('Hour of Day')
plt.ylabel('Average Total Sales')
plt.xticks(range(0, 24, 2))

plt.tight_layout()
plt.show()

# Ensure the index is datetime and properly sorted
monthly_agg.index = pd.to_datetime(monthly_agg.index)
monthly_agg = monthly_agg.sort_index()

# Loop through each drug category
for col in monthly_agg.columns:
    try:
        # Additive Decomposition
        result_add = seasonal_decompose(monthly_agg[col], model='additive', period=12)
        # Multiplicative Decomposition
        result_mul = seasonal_decompose(monthly_agg[col], model='multiplicative', period=12)

        # Plot Additive
        plt.figure(figsize=(12, 8))
        result_add.plot()
        plt.suptitle(f'Additive Decomposition of {col}', fontsize=16)
        plt.tight_layout()
        plt.show()

        # Plot Multiplicative
        plt.figure(figsize=(12, 8))
        result_mul.plot()
        plt.suptitle(f'Multiplicative Decomposition of {col}', fontsize=16)
        plt.tight_layout()
        plt.show()

    except Exception as e:
        print(f"Could not decompose {col}: {e}")

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

for col in monthly_agg.columns:
    try:
        result_add = seasonal_decompose(monthly_agg[col], model='additive', period=12)
        result_mul = seasonal_decompose(monthly_agg[col], model='multiplicative', period=12)

        # Mean and Variance of Residuals
        add_mean, add_var = result_add.resid.mean(), result_add.resid.var()
        mul_mean, mul_var = result_mul.resid.mean(), result_mul.resid.var()

        print(f"\n--- {col} ---")
        print(f"Additive Residuals: Mean = {add_mean:.2f}, Variance = {add_var:.2f}")
        print(f"Multiplicative Residuals: Mean = {mul_mean:.2f}, Variance = {mul_var:.2f}")

        # Decision Rule
        if abs(add_mean) < 0.1 and add_var < mul_var:
            print("Suggested Model: ADDITIVE (better-centered residuals, lower variance)")
        else:
            print("Suggested Model: MULTIPLICATIVE (residuals scale with trend)")
    except Exception as e:
        print(f"Error with {col}: {e}")

--- M01AB ---
Additive Residuals: Mean = -0.38, Variance = 147.70
Multiplicative Residuals: Mean = 1.00, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- M01AE ---
Additive Residuals: Mean = 0.38, Variance = 162.95
Multiplicative Residuals: Mean = 1.00, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N02BA ---
Additive Residuals: Mean = -0.45, Variance = 131.79
Multiplicative Residuals: Mean = 0.99, Variance = 0.01
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N02BE ---
Additive Residuals: Mean = 0.29, Variance = 18716.85
Multiplicative Residuals: Mean = 1.00, Variance = 0.02
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N05B ---
Additive Residuals: Mean = 3.34, Variance = 2735.99
Multiplicative Residuals: Mean = 1.01, Variance = 0.03
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- N05C ---
Additive Residuals: Mean = -0.14, Variance = 36.64
Multiplicative Residuals: Mean = 0.99, Variance = 0.12
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- R03 ---
Additive Residuals: Mean = 0.95, Variance = 1278.52
Multiplicative Residuals: Mean = 1.00, Variance = 0.04
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

--- R06 ---
Additive Residuals: Mean = 1.21, Variance = 273.80
Multiplicative Residuals: Mean = 1.01, Variance = 0.03
Suggested Model: MULTIPLICATIVE (residuals scale with trend)

# Apply log transformation to all columns
monthly_agg_log = np.log(monthly_agg)

for col in monthly_agg_log.columns:
    try:
        result = seasonal_decompose(monthly_agg_log[col], model='additive', period=12)

        plt.figure(figsize=(12, 8))
        result.plot()
        plt.suptitle(f'Log-Transformed Seasonal Decomposition of {col}', fontsize=16)
        
        # Plot residuals as dots
        plt.subplot(414)
        plt.plot(result.resid, 'o', markersize=3, alpha=0.7)
        plt.title("Residuals (Dots)")
        plt.tight_layout()
        plt.show()
        
    except Exception as e:
        print(f"Could not decompose {col}: {e}")

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

<Figure size 1200x800 with 0 Axes>

# -------------------------
# 1. Lagged Values & ACF/PACF Analysis
# -------------------------
def analyze_acf_pacf(time_series, lags=30, title=""):
    plt.figure(figsize=(12, 5))
    plt.subplot(121)
    plot_acf(time_series, ax=plt.gca(), lags=lags)
    plt.title(f'ACF - {title}')
    plt.ylabel('Autocorrelation')
    plt.subplot(122)
    plot_pacf(time_series, ax=plt.gca(), lags=lags)
    plt.title(f'PACF - {title}')
    plt.ylabel('Partial Autocorrelation')
    plt.tight_layout()
    plt.show()

    print(f"\n### ACF and PACF Analysis - {title}")
    print("General Interpretation Guide:")
    print("- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.")
    print("- PACF: Cutoff after a few lags suggests AR component.")

# -------------------------
# 2. Log Transform Utility
# -------------------------
def safe_log_transform(series, epsilon=1e-6):
    if (series <= 0).any():
        print("⚠️ Time series contains non-positive values. Adding a small constant before log transformation.")
        return np.log(series + epsilon)
    return np.log(series)

def reverse_safe_log_transform(log_series, epsilon=1e-6, was_adjusted=False):
    if was_adjusted:
        return np.exp(log_series) - epsilon
    else:
        return np.exp(log_series)

# -------------------------
# 3. ARIMA Order Selection
# -------------------------
def select_arima_order(pacf_lags, acf_lags, is_stationary=False):
    p = len(pacf_lags)
    q = len(acf_lags)
    d = 0 if is_stationary else 1
    print("\n### ARIMA Order Selection")
    print(f"Based on ACF/PACF, suggested ARIMA orders: p={p}, d={d}, q={q}")
    return p, d, q

# -------------------------
# 4. Fit ARIMA Model
# -------------------------
def fit_arima_model(time_series, order=(1, 1, 1), apply_log=True):
    used_log = False
    was_adjusted = False
    transformed_series = time_series
    model_fit = None
    lb_q_stat = np.nan
    lb_p_value = np.nan

    if apply_log:
        try:
            transformed_series = safe_log_transform(time_series)
            used_log = True
            if (time_series <= 0).any():
                was_adjusted = True
        except ValueError as e:
            print(f"⚠️ Error during log transformation: {e}")
            print("⚠️ Proceeding without log transformation.")

    try:
        model = ARIMA(transformed_series, order=order)
        model_fit = model.fit()
        print("\n### ARIMA Model Summary")
        print(model_fit.summary())

        # Augmented Dickey-Fuller Test (Stationarity of Residuals)
        result = adfuller(model_fit.resid)
        print('\nAugmented Dickey-Fuller Test:')
        print(f'ADF Statistic: {result[0]:.3f}')
        print(f'p-value: {result[1]:.3f}')
        print('Critical Values:')
        for key, value in result[4].items():
            print(f'  {key}: {value:.3f}')

        # Normality of Residuals (Jarque-Bera Test)
        jb_test = sm.stats.jarque_bera(model_fit.resid)
        print('\nJarque-Bera Test (Normality of Residuals):')
        print(f'JB Statistic: {jb_test[0]:.3f}')
        print(f'p-value: {jb_test[1]:.3f}')

        # Autocorrelation of Residuals (Ljung-Box Test)
        try:
            lb_test = acorr_ljungbox(model_fit.resid, lags=[10], return_df=True)
            print('\nLjung-Box Test (Autocorrelation of Residuals):')
            
            if isinstance(lb_test, pd.DataFrame):
                if 'lb_stat' in lb_test.columns:
                    lb_q_stat = lb_test['lb_stat'].iloc[0]
                    lb_p_value = lb_test['lb_pvalue'].iloc[0]
                elif 'lb_Q' in lb_test.columns:
                    lb_q_stat = lb_test['lb_Q'].iloc[0]
                    lb_p_value = lb_test['lb_pvalue'].iloc[0]
            else:  # Very old versions return a tuple
                lb_q_stat = lb_test[0][0]
                lb_p_value = lb_test[1][0]
            
            print(f'  Q Statistic: {lb_q_stat:.3f}')
            print(f'  p-value: {lb_p_value:.3f}')
            
        except Exception as e:
            print(f"⚠️ Error during Ljung-Box test: {e}")
            print("  Ljung-Box test results not available.")
            lb_q_stat = np.nan
            lb_p_value = np.nan

        return model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value

    except Exception as e:
        print(f"⚠️ Error fitting ARIMA model: {e}")
        return None, False, False, np.nan, np.nan

# -------------------------
# 5. Compare ARIMA with Exp. Smoothing & Prophet
# -------------------------
def compare_models(hourly_df, arima_model_fits, test_size=0.2, smoothing_level=0.2):
    print("\n### Model Comparison for All Pharmaceutical Columns")
    results = {}
    for drug in ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']:
        if drug in hourly_df:
            time_series = hourly_df[drug].dropna()
            if not time_series.empty:
                print(f"\n--- Comparison for {drug} ---")
                train = time_series[:int(len(time_series)*(1 - test_size))]
                test = time_series[int(len(time_series)*(1 - test_size)):]

                # SES
                model_ses = SimpleExpSmoothing(train).fit(smoothing_level=smoothing_level)
                forecast_ses = model_ses.forecast(steps=len(test))

                # ARIMA
                arima_mae = np.nan  # Initialize
                arima_lower_conf = np.nan
                arima_upper_conf = np.nan
                arima_lb_q_stat = np.nan
                arima_lb_p_value = np.nan

                if drug in arima_model_fits and arima_model_fits[drug][0] is not None:
                    model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value = arima_model_fits[drug]  # Get Ljung-Box results
                    arima_forecast = model_fit.get_forecast(steps=len(test)).predicted_mean
                    conf_int = model_fit.get_forecast(steps=len(test)).conf_int(alpha=0.05)
                    arima_lower_conf = conf_int.iloc[:, 0].values
                    arima_upper_conf = conf_int.iloc[:, 1].values
                    arima_lb_q_stat = lb_q_stat
                    arima_lb_p_value = lb_p_value

                    if used_log:
                        arima_forecast = reverse_safe_log_transform(arima_forecast, was_adjusted=was_adjusted)
                        if not np.isnan(arima_lower_conf).any() and not np.isnan(arima_upper_conf).any():
                            arima_lower_conf = reverse_safe_log_transform(arima_lower_conf, was_adjusted=was_adjusted)
                            arima_upper_conf = reverse_safe_log_transform(arima_upper_conf, was_adjusted=was_adjusted)
                        else:
                            arima_lower_conf = np.nan
                            arima_upper_conf = np.nan

                    # ARIMA Evaluation
                    y_true_arima = np.array(test)
                    y_pred_arima = arima_forecast[:len(test)]  # Ensure lengths match
                    if len(y_true_arima) > 0 and len(y_pred_arima) > 0:
                        arima_mae = mean_absolute_error(y_true_arima, y_pred_arima)

                    # ARIMA Confidence Intervals
                    print("\nARIMA Confidence Intervals:")
                    if not np.isnan(arima_lower_conf).any() and not np.isnan(arima_upper_conf).any():
                        for i in range(min(5, len(test))):  # Print a few examples
                            print(f"  Time {test.index[i]}: Forecast={arima_forecast[i]:.4f}, Lower={arima_lower_conf[i]:.4f}, Upper={arima_upper_conf[i]:.4f}")
                    else:
                        print("  Confidence intervals not available (NaNs present).")
                    
                    # ARIMA Ljung-Box Test Results
                    print("\nARIMA Ljung-Box Test Results:")
                    print(f"  Q-Statistic: {arima_lb_q_stat:.4f}")
                    print(f"  P-Value: {arima_lb_p_value:.4f}")

                else:
                    print(f"⚠️ ARIMA model not found or failed to fit for {drug}. Skipping ARIMA evaluation.")
                    arima_forecast = np.full(len(test), np.nan)

                # Prophet
                try:
                    prophet_df = pd.DataFrame({'ds': train.index, 'y': train.values})
                    model_prophet = Prophet(daily_seasonality=True)
                    model_prophet.fit(prophet_df)
                    future = pd.DataFrame({'ds': test.index})
                    forecast_prophet = model_prophet.predict(future)
                    prophet_mean = forecast_prophet['yhat'].values
                    prophet_lower = forecast_prophet['yhat_lower'].values
                    prophet_upper = forecast_prophet['yhat_upper'].values
                    prophet_mae = mean_absolute_error(test, prophet_mean)

                    # Prophet Prediction Intervals
                    print("\nProphet Prediction Intervals:")
                    for i in range(min(5, len(test))):  # Print a few examples
                        print(f"  Time {test.index[i]}: Forecast={prophet_mean[i]:.4f}, Lower={prophet_lower[i]:.4f}, Upper={prophet_upper[i]:.4f}")

                except Exception as e:
                    print(f"⚠️ Error fitting or predicting with Prophet for {drug}: {e}")
                    prophet_mean = np.full(len(test), np.nan)
                    prophet_lower = np.full(len(test), np.nan)
                    prophet_upper = np.full(len(test), np.nan)
                    prophet_mae = np.nan

                # SES Evaluation
                ses_mae = mean_absolute_error(test, forecast_ses)
                
                results[drug] = {
                    "MAE_ARIMA": arima_mae,
                    "MAE_ExpSmoothing": ses_mae,
                    "MAE_Prophet": prophet_mae,
                    "ARIMA_LjungBox_Q": arima_lb_q_stat,
                }
                
                best_model = None
                model_maes = [("ARIMA", arima_mae), ("Exponential Smoothing", ses_mae), ("Prophet", prophet_mae)]
                valid_models = [(name, mae) for name, mae in model_maes if not np.isnan(mae)]
                
                print(f"  - MAE ARIMA: {arima_mae:.4f}")
                print(f"  - MAE Exponential Smoothing: {ses_mae:.4f}")  # Fixed variable name here
                print(f"  - MAE Prophet: {prophet_mae:.4f}")
                
                if valid_models:
                    best_model = min(valid_models, key=lambda item: item[1])[0]
                    print(f"  - Best Model: {best_model}")
                else:
                    print("  - All models failed to evaluate.")

            else:
                print(f"  - No data to process for {drug}")
        else:
            print(f"  - Column '{drug}' not found in DataFrame.")
    return results

# -------------------------
# 6. Pipeline Execution
# -------------------------
arima_model_fits = {}

for drug in ['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']:
    if drug in hourly_df:
        time_series = hourly_df[drug].dropna()
        print(f"\n### ACF/PACF Analysis for {drug}")
        # Analyze ACF/PACF on the original data first
        analyze_acf_pacf(time_series, lags=30, title=f'{drug} (original)')

        # Then analyze on the potentially log-transformed data
        analyzed_series = safe_log_transform(time_series)
        analyze_acf_pacf(analyzed_series, lags=30, title=f'{drug} (log-transformed)')

        pacf_significant_lags = [1]
        acf_significant_lags = []
        is_stationary = False

        p, d, q = select_arima_order(pacf_significant_lags, acf_significant_lags, is_stationary)
        model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value = fit_arima_model(time_series, order=(p, d, q), apply_log=True)
        if model_fit is not None:  # Only store if model fit successfully
            arima_model_fits[drug] = (model_fit, used_log, was_adjusted, lb_q_stat, lb_p_value)
        else:
            arima_model_fits[drug] = (None, used_log, was_adjusted, np.nan, np.nan)  # Store None to indicate failure
    else:
        print(f"Column '{drug}' not found in DataFrame.")

if arima_model_fits:
    comparison_results = compare_models(hourly_df, arima_model_fits)
    print("\nFinal Comparison Results:")
    results_df = pd.DataFrame.from_dict(comparison_results, orient='index')
    results_df = results_df.rename(columns={'MAE_ARIMA': 'MAE ARIMA', 'MAE_ExpSmoothing': 'MAE Exp. Smoothing', 'MAE_Prophet': 'MAE Prophet',
                                      'ARIMA_LjungBox_Q': 'ARIMA Ljung-Box Q', 'ARIMA_LjungBox_pvalue': 'ARIMA Ljung-Box p-value'})
    # Determine the best model for each drug
    def get_best_model(row):
        row_filtered = row[row.notna()]
        if not row_filtered.empty:
            return row_filtered.idxmin().replace('MAE ', '')
        return "N/A"
    results_df['Best Model'] = results_df[['MAE ARIMA', 'MAE Exp. Smoothing', 'MAE Prophet']].apply(get_best_model, axis=1)
    print(results_df.to_markdown())
else:
    print("No ARIMA models were fit, so comparison was skipped.")

### ACF/PACF Analysis for M01AB

### ACF and PACF Analysis - M01AB (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - M01AB (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  M01AB   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -161937.782
Date:                Wed, 07 May 2025   AIC                         323879.564
Time:                        05:51:52   BIC                         323897.225
Sample:                    01-02-2014   HQIC                        323885.093
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4801      0.003   -143.032      0.000      -0.487      -0.474
sigma2        35.5727      0.180    197.822      0.000      35.220      35.925
===================================================================================
Ljung-Box (L1) (Q):                1030.89   Jarque-Bera (JB):              4522.69
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.07   Skew:                             0.48
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.10
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.848
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 4520.699
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5479.990
  p-value: 0.000

### ACF/PACF Analysis for M01AE

### ACF and PACF Analysis - M01AE (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - M01AE (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  M01AE   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -161162.018
Date:                Wed, 07 May 2025   AIC                         322328.035
Time:                        05:52:10   BIC                         322345.696
Sample:                    01-02-2014   HQIC                        322333.565
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4762      0.003   -139.780      0.000      -0.483      -0.469
sigma2        34.4948      0.177    194.904      0.000      34.148      34.842
===================================================================================
Ljung-Box (L1) (Q):                 925.40   Jarque-Bera (JB):              3777.83
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.06   Skew:                             0.44
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.01
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -49.021
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 3777.922
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5108.978
  p-value: 0.000

### ACF/PACF Analysis for N02BA

### ACF and PACF Analysis - N02BA (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - N02BA (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  N02BA   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -158341.704
Date:                Wed, 07 May 2025   AIC                         316687.409
Time:                        05:52:25   BIC                         316705.070
Sample:                    01-02-2014   HQIC                        316692.938
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.4774      0.003   -151.476      0.000      -0.484      -0.471
sigma2        30.8534      0.143    215.770      0.000      30.573      31.134
===================================================================================
Ljung-Box (L1) (Q):                1018.04   Jarque-Bera (JB):              9182.37
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               0.84   Skew:                             0.61
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.69
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.343
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 9183.360
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 5499.136
  p-value: 0.000

### ACF/PACF Analysis for N02BE

### ACF and PACF Analysis - N02BE (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ACF and PACF Analysis - N02BE (log-transformed)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.

### ARIMA Order Selection
Based on ACF/PACF, suggested ARIMA orders: p=1, d=1, q=0
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

### ARIMA Model Summary
                               SARIMAX Results                                
==============================================================================
Dep. Variable:                  N02BE   No. Observations:                50532
Model:                 ARIMA(1, 1, 0)   Log Likelihood             -167559.956
Date:                Wed, 07 May 2025   AIC                         335123.913
Time:                        05:52:40   BIC                         335141.574
Sample:                    01-02-2014   HQIC                        335129.442
                         - 10-08-2019                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.3998      0.004   -104.461      0.000      -0.407      -0.392
sigma2        44.4353      0.239    186.186      0.000      43.968      44.903
===================================================================================
Ljung-Box (L1) (Q):                 274.49   Jarque-Bera (JB):              1559.15
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.03   Skew:                            -0.16
Prob(H) (two-sided):                  0.03   Kurtosis:                         3.80
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Augmented Dickey-Fuller Test:
ADF Statistic: -48.837
p-value: 0.000
Critical Values:
  1%: -3.430
  5%: -2.862
  10%: -2.567

Jarque-Bera Test (Normality of Residuals):
JB Statistic: 1559.593
p-value: 0.000

Ljung-Box Test (Autocorrelation of Residuals):
  Q Statistic: 2389.539
  p-value: 0.000

### ACF/PACF Analysis for N05B

### ACF and PACF Analysis - N05B (original)
General Interpretation Guide:
- ACF: Slow decay suggests trend/non-stationarity; cutoff after a few lags suggests MA component.
- PACF: Cutoff after a few lags suggests AR component.
⚠️ Time series contains non-positive values. Adding a small constant before log transformation.

	datum	M01AB	M01AE	N02BA	N02BE	N05B	N05C	R03	R06
0	1/5/2014	14.00	11.67	21.3	185.95	41.0	0.0	32.0	7.0
1	1/12/2014	29.33	12.68	37.9	190.70	88.0	5.0	21.0	7.2
2	1/19/2014	30.67	26.34	45.9	218.40	80.0	8.0	29.0	12.0
3	1/26/2014	34.00	32.37	31.5	179.60	80.0	8.0	23.0	10.0
4	2/2/2014	31.02	23.35	20.7	159.88	84.0	12.0	29.0	12.0

Forecasting Point-of-Sale Trends in Pharmaceuticals Using ARIMA¶

Sasha Nabila Fortuna¶

Background Problem¶

Objectives¶

Dataset Description¶

Methodology¶

Import Libraries¶

Load Dataset¶

Exploratory Data Analysis¶

Comparison Drug Products Over Time¶

Average Sales by Day for All Products¶

Average Hourly Sales for All Products¶

Time Series Decomposition¶

ARIMA Modelling¶

Conclusion¶

Further Improvement¶

References¶

	datum	M01AE	N02BA	N02BE	N05B	R06	Year	Month	Hour	Weekday Name
0	1/2/2014 8:00	0.67	0.4	2.0	0.0	1.0	2014	1	8	Thursday
1	1/2/2014 9:00	0.00	1.0	0.0	2.0	0.0	2014	1	9	Thursday
2	1/2/2014 10:00	0.00	0.0	3.0	2.0	0.0	2014	1	10	Thursday
3	1/2/2014 11:00	0.00	0.0	2.0	1.0	0.0	2014	1	11	Thursday
4	1/2/2014 12:00	2.00	0.0	5.0	2.0	0.0	2014	1	12	Thursday

	datum	M01AB	M01AE	N02BA	N02BE	N05B	N05C	R03	R06	Year	Month	Hour	Weekday Name
0	1/2/2014	0.0	3.67	3.4	32.40	7.0	0.0	0.0	2.0	2014	1	248	Thursday
1	1/3/2014	8.0	4.00	4.4	50.60	16.0	0.0	20.0	4.0	2014	1	276	Friday
2	1/4/2014	2.0	1.00	6.5	61.85	10.0	0.0	9.0	1.0	2014	1	276	Saturday
3	1/5/2014	4.0	3.00	7.0	41.10	8.0	0.0	3.0	0.0	2014	1	276	Sunday
4	1/6/2014	5.0	1.00	4.5	21.70	16.0	2.0	6.0	2.0	2014	1	276	Monday